Global gene expression profiling under nitrogen stress identifies key genes involved in nitrogen stress adaptation in maize (Zea mays L.)

Maize is a heavy consumer of fertilizer nitrogen (N) which not only results in the high cost of cultivation but may also lead to environmental pollution. Therefore, there is a need to develop N-use efficient genotypes, a prerequisite for which is a greater understanding of N-deficiency stress adaptation. In this study, comparative transcriptome analysis was performed using leaf and root tissues from contrasting inbred lines, viz., DMI 56 (tolerant to N stress) and DMI 81 (susceptible to N stress) to delineate the differentially expressed genes (DEGs) under low-N stress. The contrasting lines were grown hydroponically in modified Hoagland solution having either sufficient- or deficient-N, followed by high-throughput RNA-sequencing. A total of 8 sequencing libraries were prepared and 88–97% of the sequenced raw reads were mapped to the reference B73 maize genome. Genes with a p value ≤ 0.05 and fold change of ≥ 2.0 or ≤ − 2 were considered as DEGs in various combinations performed between susceptible and tolerant genotypes. DEGs were further classified into different functional categories and pathways according to their putative functions. Gene Ontology based annotation of these DEGs identified three different functional categories: biological processes, molecular function, and cellular component. The KEGG and Mapman based analysis revealed that most of the DEGs fall into various metabolic pathways, biosynthesis of secondary metabolites, signal transduction, amino acid metabolism, N-assimilation and metabolism, and starch metabolism. Some of the key genes involved in N uptake (high-affinity nitrate transporter 2.2 and 2.5), N assimilation and metabolism (glutamine synthetase, asparagine synthetase), redox homeostasis (SOD, POX), and transcription factors (MYB36, AP2-EREBP) were found to be highly expressed in the tolerant genotype compared to susceptible one. The candidate genes identified in the present study might be playing a pivotal role in low-N stress adaptation in maize and hence could be useful in augmenting further research on N metabolism and development of N-deficiency tolerant maize cultivars.

Nitrogen (N) is the most important macronutrient required for the growth and development of crop plants. It is a crucial structural element in major biomolecules, viz., chlorophyll, proteins, enzymes, nucleic acids, and various hormones. Besides nutrients, it also acts as a signal molecule. In the soil, it is available in the more common water-soluble nitrate (NO 3 − ) form, relatively less common ammonium (NH 4 + ) form and to a lesser extent as proteins, peptides, or amino acids 1 . Being an essential nutrient, N-deficiency causes chlorosis in leaves, especially in lower leaves, restricts bud growth, and reduces overall plant growth [2][3][4] . Its deficiency at the early vegetative stage of plant life adversely affects crop yield which cannot be reversed by applying N at later stages 5 . Since its availability strongly affects crop productivity, a vast amount of N fertilizers is applied to increase crop yield. However, nitrogen use efficiency (NUE) in cereal crops ranges from 40 to 50% 6,7 . Thus, a significant amount of fertilizer N is leached into the groundwater or evaporated into the environment and thereby, contributes to contamination of ground and surface water, emission of greenhouse gase i.e. nitrous oxide, and also deteriorates Scientific Reports | (2022) 12:4211 | https://doi.org/10.1038/s41598-022-07709-z www.nature.com/scientificreports/ soil health [8][9][10][11][12] . Hence, developing crop genotypes with improved NUE would be of prime importance to achieve sustainable agriculture and high productivity under low input conditions with low environmental footprint. NUE involves efficiency in N uptake, assimilation, remobilization, and utilization. N uptake and metabolism pathways in plants have been well elucidated. The nitrate transporters present in plant root cells-high (NRT1) and low (NRT2) affinity transporters-help in N uptake from the soil, which is then further metabolized to nitrite and ammonium by nitrate reductase (NR) and nitrite reductase (NiR), respectively. This converted ammonium is then incorporated into the organic form (amino acids) by glutamine synthetase (GS) and glutamate synthase/glutamine-2-oxoglutarate aminotransferase (GOGAT) enzymes, also known as GS/GOGAT cycle 13 . In higher plants, there are two types of glutamine synthetase: GS1 isoenzyme, cytosolic form (have 5 isoforms in Arabidopsis), and GS2 isoenzyme, present in the plastid. Similarly, two types of glutamate synthase are also reported-ferredoxin-dependent glutamate synthase (Fd-GOGAT), present in the chloroplast in shoots, and NADH-GOGAT, present in root plastids 14 . The GS2 and Fd-GOGAT assimilate ammonium into glutamine and further into glutamate, respectively in the chloroplast, and are essential for survival under photorespiratory conditions 15 . The GS1 assimilates ammonium into glutamine in root cytosol 16,17 . Assimilated N is transported as asparagine, glutamine, aspartate, and glutamate for storage, assimilation, and utilization 18 . The plants' ability to effectively remobilize N into maturing fruits or grains is very important to NUE 19 , especially in cereal crops where the grain is also economically important.
NUE is a complex trait and to date considerable efforts have been made to understand the molecular basis of plant responses to N and identifying N-responsive genes, and regulatory factors so that their expression could be manipulated for better NUE [20][21][22][23] . For instance, in Arabidopsis, microarray studies revealed that expression levels of various genes vary with different concentrations of nitrate both for long-term, and short-term basis 24,25 . In rice, expression analysis of 10,422 genes by microarray revealed a significant difference in the expression level of 471 genes in root tissue 26 . Similarly, in maize, a few studies have been carried out to identify N stress-responsive genes [27][28][29][30] . However, a major limitation in most of these studies was that these have been carried out by using a single genotype. Without comparing the transcriptional differences between N-stress tolerant and susceptible genotypes, it is not prudent to separate N-stress tolerant genes from stress-responsive genes.
Maize (Zea mays L.) is an important cereal crop for feed, food and industrial raw material. It is the most produced grain in the world. Predominantly, hybrid maize cultivars are cultivated in intensive cropping systems, with high external N input. It is a heavy N consumer and it is highly susceptible to N stress especially in the vegetative stage when uptake and utilization of N are at their peak. To develop N efficient maize genotypes, it is highly essential to delineate the candidate genes and master regulators playing a critical role in NUE in maize. To the best of our knowledge, there are two reports in which contrasting genotypes were used for identification of NUE genes: (1) Chen et al. 31 , has reported gene expression changes in response to low nitrogen stress in leaf tissues of contrasting maize inbred lines using the Affymetrix maize genome array, (2) Zamboni et al. 32 , has analyzed transcriptional changes in the root of a high and a low nitrogen use efficient maize inbred line in the response to nitrate treatment for 24 h in 7-day old seedling. Both of these studies have their limitations, viz., in the first one, differentially expressed genes (DEGs) in leaf tissue was studied, but not in the root tissue while in the second one, expression changes in root tissue of 7-day old seedling was studied and with a short duration of nitrate treatment. However, understanding of the NUE trait can be further improved by identifying genes in leaf as well as in root tissues at a time under effective treatment (longer duration) for nitrogen stress. Hence, the present study aimed to identify key genes involved in determining nitrogen use efficiency in maize. For this, gene expression was analyzed in root and leaf tissue from contrasting (N-deficiency tolerant and susceptible) tropical maize inbred lines under N-deficiency conditions vis-a-vis control conditions using high-throughput RNA sequencing (RNA-seq) technology.

Results
Phenotypic performance of maize genotypes. Under N-deficiency conditions, the shoot biomass decreased by 56.3% and 68.2% in the N-stress tolerant (DMI 56) and susceptible genotype (DMI 81), respectively, while the root biomass increased significantly in both the genotypes (Table 1, Fig. 1, Supplementary Fig  S1). The root of DMI 56 was longer (by 110.8%) when it was grown under low-N for 21 days. The susceptible genotype, though having initially longer roots as compared to the tolerant genotype, could not sufficiently expand its root under N-deficiency (only a 24% increase in root length). Under N-sufficient conditions, the root system of the tolerant genotype was less extensive than those of the susceptible genotype although there was no major difference in the shoot. However, under N-deficiency conditions, the root length of the tolerant genotype was appreciably higher (by 26

Identification of differentially expressed genes under nitrogen-deficiency.
To identify significant DEGs under N-deficiency stress in tropical maize, transcript profiling from contrasting inbred lines was studied under N-deficient and control conditions (see ***"Methods" section for details). As the correct understanding of the differential expression of genes is key to infer phenotypic variations observed among genotypes, therefore, 6 comparisons were made between DMI 56 (tolerant line) and DMI 81(susceptible line) root and leaf tissues to analyze DEGs. A total of 1908, 2444, 1827, 2521, 1873, and 1034 significant DEGs were mapped and annotated in the case of 1st, 2nd 3rd, 4th, 5th, and 6th combinations, respectively ( Table 2). A complete list of significant DEGs is provided in Supplementary Table S1 (excel file). The chromosome-wise distribution of DEGs from various combinations was visualized using Circos representation for all the ten chromosomes (Fig. 2). A total of seven and two DEGs in root were common in all three down-regulated and up-regulated combinations, respectively (Fig. 3). Similarly, in leaf, five DEGs were common in all three down-regulated combinations whereas one DEG was found common in up-regulated combinations.
The MapMan based visualization of the expression of DEGs onto metabolic pathways revealed that the maximum number of up-regulated genes were related to secondary metabolism followed by cell wall and lipids processes in DMI 81 leaf as compared to DMI 56 leaf under low-N stress (Fig. 4A). Other than C-1 metabolism, most of DEGs related to photosynthesis, starch-sucrose metabolism, N metabolism were down-regulated. While in the case of root under low-N stress, the maximum number of DEGsbelonged to the cell wall, lipids, and secondary metabolism pathways and most of them showed down-regulation. Besides, DEGs related to NO 3 metabolism, amino-acid metabolism, and photosynthesis were down-regulated and C-1 metabolism was upregulated (Fig. 4B). Further, KEGG Pathway analysis revealed that the maximum number of the DEGs could be classified into metabolic pathways, biosynthesis of secondary metabolites, phenylpropanoid biosynthesis in both the comparisons i.e. DMI 81 leaf and root as compared to DMI 56 leaf and root, respectively, under low-N stress (Fig. 5).
WEGO plots represent up-regulated and down-regulated GO annotations which were distributed into three categories namely molecular functions, cellular components, and biological processes. Under low-N stress, a   Table S1 given as excel file). Among the down-regulated genes, most of the genes belonged to biological processes and cellular components followed by molecular function (Fig. 6A). A similar pattern was observed in up-regulated genes (Fig. 6B). When we analyzed the DEGs in the root under low-N stress, 1056 genes were found up-regulated, whereas 1388 were down-regulated in DMI 81 as compared to DMI 56 (Supplementary Table S1). Genes involved in the cellular category and followed by biological process were prominently differentially expressed in the root tissue (Fig. 7). Further, top 20 up-regulated and down-regulated genes under N-deficiency in root and leaf tissues of susceptible genotype in comparison to tolerant genotype were shortlisted (Figs. 8,9). In all combinations, few of the DEGs were uncharacterized which means that their sequence does not match any annotated genes in the database, hence they can be called novel transcripts. Furthermore, the significant DEGs in leaf and root of DMI   (Fig. 10). These key DEGs selected for validation encode genes and transcription factors playing a pivotal role in nitrogen metabolism, transport, and signaling mechanisms. For example, the asparagine synthetase enzyme helps in ammonium assimilation and also plays an important role in nitrogen assimilation, recycling, transport, and storage in plants. Amongst the other DEGs, HAT 2.3 encodes for a high-affinity transporter for nitrogen and nitrate transporter 2.5 is involved in the constitutive high-affinity transport system under long-term N starvation conditions in plants.

Discussion
The present study was undertaken with the aim to identify key genes playing crucial roles in N stress tolerance in tropical maize. For this, we studied the transcriptome of leaf and root tissues from contrasting maize genotypes under N-deficiency as well as control conditions. Maize plants grown hydroponically under low-N stress conditions exhibited visual symptoms of N-deficiency such as stunted growth, pinkish-red coloration in shoots, yellow coloration in older leaves, upright leaves with light green/yellow color, and burnt leaf margin, etc. These symptoms were more prominent in the susceptible genotype as compared to the tolerant one ( Supplementary  Fig S1). Soltabayeva et al. 33 have reported accelerated yellowing and senescence of old leaves as one of the typical  (Table 1). This indicates stress-related growth retardation, highlighting the prominent role of N for biomass accumulation. However, there was a pronounced increase in root length under low-N stress conditions as compared to N-sufficient conditions ( Fig. 1; Table 1). Importantly, the percentage increase in root length in the tolerant genotype was much higher than the susceptible genotype. The increase in root length indicates the adaptive response of maize plants under N-limitation to maximize N uptake by increasing the surface area for acquisition. Further, it suggests that the tolerant genotype Transcriptome analysis and genes responsible for nitrogen-deficiency stress adaptation. Transcriptomics is an important and powerful approach being used for global gene expression profiling in different crop plants under abiotic, biotic, and nutrient deficiency stresses [21][22][23]35,36 . In the present study, transcriptomic comparison of low-N stress-tolerant (DMI 56) and susceptible (DMI 81) genotypes was attempted to delineate potential candidate genes involved in N-deficiency stress in tropical maize. The cDNA samples from root and leaf tissues were sequenced and 40,724-41,868 total genes were identified, out of which significant genes ranged from 1034 to 2521 in six different comparisons ( Table 2). The number of total genes and significant genes identified in the present study were comparable to other studies 31,37,38 . The number of DEGs in both root and leaf tissues were higher in the tolerant genotype (DMI 56) as compared to the susceptible genotype (DMI 81) under N stress (Table 2). Further, it was observed that in the tolerant genotype, the number of up-regulated genes was more than the down-regulated genes as compared to the susceptible genotype (Fig. 3). In various comparisons, the number of down-regulated genes ranged from 827 to 1388 in root and 443 to 1354 in leaf while the number of up-regulated genes ranged from 507 to 1081 in root and 591 to 1167 in leaf ( Fig. 3; Supplementary  Table S1). The up-regulated and down-regulated genes obtained in our study are comparable to those reported recently in Tibetan hulless barley 39 and sorghum 40 but lower than rice 41 under N stress.
Our study revealed that most of the DEGs were mainly confined to secondary metabolism, cell-wall, and lipid component in all six combinations ( Fig. 4; Supplementary Figs. S2, S3). Some DEGs also showed significant up-regulation and down-regulation in C-1 metabolism, NO 3 metabolism, starch sucrose cycle, photosynthesis, amino-acid metabolism, and photorespiration. Interestingly, most of the DEGs were up-regulated and downregulated in root tissues of DMI 56 and DMI 81, respectively, under low-N stress as compared to their respective control ( Supplementary Figs. S2, S3). Further, the number of DEGs was more in root than in the leaf. In Arabidopsis thaliana, DEGs associated with the above-mentioned processes played a very important role in adaptation   42 . Similarly, our finding also corroborates with previously identified DEGs related to various pathways in rice 43 and wheat under varied nitrogen supplies 44 . Different studies have reported the involvement of genes related to amino acid metabolism, lipid metabolism, energy metabolism, and signal transduction pathways for tolerance against N stress [45][46][47] . For example, in sorghum roots, under low-N stress, DEGs associated with amino-acid metabolism play a key role 40 . Similarly, protein kinases (PK) are prominently involved in response to N stress in Arabidopsis thaliana leaves and roots 44 . In this study, genes related to metabolic pathways, biosynthesis of secondary metabolite, plant hormone signal, biosynthesis of amino acids, photosynthesis, MAPK signaling, etc. were identified in both root and leaves ( Fig. 5; Supplementary Figs. S4, S5). The present findings support the results from previous studies and confirm that all these pathways form a complex network in N-stress adaptation in plants.    54 and Arabidopsis under heat stress 55 . CDP-diacylglycerol-serine O-phosphatidyltransferase catalyzes base-exchange reaction where phosphatidylcholine is replaced by serine and is involved in amino acid metabolism 56 . Some of the prominent down-regulated genes in susceptible line root compared to tolerant line root are Putative cytochrome P450 protein (Zm00001d053586), protein ALP1-like isoform X1 (LOC103645936), Serine/arginine-rich splicing factor RS2Z32 (Zm00001d037543), Katanin p80 WD40 repeat-containing subunit B1 homolog (Zm00001d032598), Glycosyltransferases (Zm00001d038981), Peptidylprolyl isomerase ( Table 3). The cytochrome P450 (CYP) is involved in various metabolic pathways and performs an important function in metabolic reactions leading to accumulation of secondary metabolites that protect the plant from many biotic www.nature.com/scientificreports/ and abiotic stresses 41,57 . Thus, down-regulation of this gene under low-N stress may contribute to the susceptible nature of any genotype while up-regulation may lead to tolerant nature. Protein ALP1-like isoform X1 is a type of transposase having nuclease activity and has been reported to provide tolerance under drought stress 58 . SR (serine/arginine-rich) proteins are a highly conserved family of RNA-binding proteins which play a crucial role in pre-mRNA splicing. This RNA splicing has been shown to be tissue-specific, developmentally regulated, and stress-responsive 59 . Katanin is a microtubule-severing protein that regulates cell division and its orientation during plant growth. Its working is under hormonal control 60 . Glycosyltransferase has been demonstrated to help in maintaining membrane integrity under abiotic stress conditions especially during chilling stress 61 . Peptidylprolyl isomerase is a domain in cyclophilins, a ubiquitous protein, which is involved in a wide range of cellular processes such as signaling, cell division, transcriptional regulation under various abiotic stresses 62 . Since these genes have been shown to involve in abiotic or biotic stress tolerance previously in different plants, therefore their up-regulation might play important role in low-N stress tolerance. Similarly, few genes showing maximum up-regulation in the susceptible line leaf tissue compared to tolerant line leaf tissue under low-N stress are the Nodulin-related protein 1 (Zm00001d008397), LRR receptor-like serine/threonine-protein kinase (Zm00001d051093), Cysteine-rich receptor-like protein kinase (Zm00001d008488),   [64][65][66] . Similarly, Cysteine-rich receptor-like protein kinases has been shown to be up-regulated under N-deficiency stress in rice 41 . In our study, transposases showed significant up-regulation in root and leaf tissues in susceptible line and down-regulation in the tolerant line under low-N stress conditions. In Arabidopsis thaliana, FIP1 (FtsH5 Interacting Protein) protein expression is modulated by light stress and it has been reported that under high light intensity, oxidative, salt, and osmotic stress its expression is downregulated 67 . In our study, the DEGs exhibiting maximum down-regulation in leaf tissue under low-nitrogen stress are Ferredoxin (Zm00001d049732), Cytochrome P450 protein, Serine carboxypeptidase-like 19 (Zm00001d003530), Transcription factor MYB36 (Zm00001d005784), Peptidylprolylisomerase (Table 4; Fig. 10). Ferredoxin is related to nitrate/nitrite assimilation, ferredoxin reduction, and the pentose phosphate pathway and has been shown to be downregulated in rice under N-deficiency stress 41 . Further, the top 20 up-and down-regulated DEGs in leaf and root tissues of DMI 56 and DMI 81 under low-N stress conditions compared to respective control are given in Supplementary Table S2-S9. Gene ontology analysis showed several DEGs involved in GO terms such as nutrient reservoir activity, molecular carrier activity, detoxification, and immune system process were significantly up-regulated in DMI 81 leaf compared to DMI 56 leaf (Fig. 6). While in DMI 81 root, DEGs involved in symplast and cell junction GO terms were down-regulated while detoxification and immune system GO terms were up-regulated Figure 10. Comparison of expression analysis of selected nitrogen stress-responsive genes via qRT-PCR (represented by blue colour) and NGS approach (represented by red colour) in maize inbreds, viz., DMI 56 (tolerant to nitrogen stress) and DMI 81 (susceptible to nitrogen stress) in response to nitrogen stress treatment. 1-6 numbers on X-axis represent comparisons in which a particular gene has significant expression. 1-6 corresponds to DMI 81_SN+ v/s DMI 81_SN−, DMI 81_RN+ v/s DMI 81_RN−, DMI 56_SN+ v/s DMI 56_SN−, DMI 56_RN+ v/s DMI 56_RN−, DMI 56_SN− v/s DMI 81_SN− and DMI 56_RN− v/s DMI 81_RN−, respectively. Letter S, R, N+, and N− represents leaf, root, sufficient nitrogen, and nitrogen-deficiency conditions respectively. Y-axis represents the log2 fold change in expression level. www.nature.com/scientificreports/ (Fig. 7). Further, genes related to amino-acid synthesis and metabolism (for example; Zm00001d048050, Zm00001d028750) were down-regulated, while genes related to the release of amino acids were up-regulated (for example; Zm00001d009944, Zm00001d031486, Zm00001d016476, and Zm00001d012667) which suggests that during the early stages of N-deficiency, recycling might be the main source of fulfilling the initial demand of N for growth and development of plant 41 . Previously, many transcription factors (TFs), such as WRKY, MYB, bHLH, bZIP have been reported to play an important role under N-deficiency stress 68,69 . In our study, few members of WRKY and bZIP TF families were found to be up-regulated in the tolerant genotype compared to susceptible one or down-regulated in the susceptible genotype compared to tolerant one. Similar results have been observed in rice where WRKY TFs were up-regulated under N-deficiency stress 23 . Besides, many other members belonging to various TFs such as basic helix loop helix, dehydrin, and late embryogenic abundant TF were found to be down-regulated in the root tissue in the current investigation. These TFs might plays important role in signaling and conferring tolerance against low-N stress in maize. Overall, it was observed that differential expression of genes was more prominent in roots as compared to leaf. This may stem from the fact that roots are generally responsible for nutrient acquisition (such as N) and under a nutrient-deficient scenario, more biological activity, like signal transduction, etc., would occur more in roots. Understanding the differential expression profiles of key genes in the roots of contrasting lines is important in elucidating the genetic determinants of N-deficiency tolerance in maize.

Conclusion
Comparative transcriptome analysis reveals that adaptive characters like increased root length and decreased shoot biomass under N-deficient conditions are also linked to the dysregulation of various genes involved in different metabolic pathways. Several potential genes were identified which might be involved in determining NUE in maize. The genes encoding for N transporters, enzymes involved in amino acid metabolism, TFs (MYB 36, AP2-EREBP, etc.), and stress-responsive-genes exhibited the genotypic dependent pattern of expression. The present study opens up the scope for the investigation to further examine and elucidate the precise role of highly dysregulated key genes in N-deficiency stress tolerance in maize. Further, the candidate genes identified in the present study may be utilized for molecular marker-assisted breeding towards the development of low-N stress-tolerant maize plants.

Methods
Plant material, stress treatment, and RNA extraction. Two maize inbred lines, DMI 56 and DMI 81 were used in this study which was developed by the ICAR-Indian Institute of Maize Research, India. These lines were selected based on a previous study that found DMI 56 as a tolerant line and DMI 81 as a susceptible line for N stress under field conditions 19 . The seeds were surface sterilized with 70% ethanol for 2 min followed by washing with sterile water 4-5 times. Subsequently, the seeds were treated with 0.1% HgCl 2 for 10 min followed by washing with sterile water four times and allowed to germinate at room temperature in wet germination paper in the dark. After 1 week, germinated seedlings of almost similar length were grown hydroponically in plastic trays for 3 days in Hoagland solution having sufficient N, after which the seedlings were carefully removed and two sets were prepared using the 10-day-old seedlings for the study. One set was transferred to modified Hoagland's hydroponic solution with sufficient N while the second set was allowed to grow under N-deficient conditions. Plants were grown until symptoms of N-deficiency were observed which was after 21 days in green-house under controlled conditions, i.e. 30/22 (± 2) °C day/night temperatures, 14/10 dark/light hours and 40-50% humidity. The nutrient medium (Hoagland solution) was changed every third day. In the nutrient medium, nitrate and ammonium-containing salt like KNO 3   www.nature.com/scientificreports/ average base content per read, GC distribution in the reads, etc. RNA-seq generated a total of 63.51 GB of highquality data from eight c-DNA libraries and clean processed reads were aligned to the reference genome i.e. maize B73 version 4 using minimap2. Mapping was carried out using minimap2 software at default parameters with Zea mays ref_B73_RefGen_v4 assembly. The quality reads showed 88.10-97.55% mapping percentage to the B73 genome (Supplementary Table 10). Annotation of reference genes was carried out using Uniprot Batch retrieval and Kyoto Encyclopedia of Genes and Genomes (KEGG). Filtered reads are mapped to the genes using bwa 0.7.5 (https:// bioweb. paste ur. fr/ packa ges/ pack@ bwa@0. 7. 5a) and SAM tools.0.1.19 (https:// sourc eforge. net/ proje cts/ samto ols/ files/ samto ols/0. 1. 19/ samto ols-0. 1. 19. tar. bz2) were used to calculate counts mapped with respect to each gene in RNA-Seq data. Later on DESeq R package (https:// www. bioco nduct or. org/ packa ges//2. 10/ bioc/ html/ DESeq. html) were used to do differential expression analysis.. The details of samples and various comparisons carried out in the present study are mentioned in Table 2. Mercator was used for annotation and Mapman 3.5.1R2 for filtering differentially expressed genes (DEGs). Significant DEGs were identified based on the statistical significance (p value ≤ 0.05) and log2 Fold change (≥ 2 for up-regulated and ≤ − 2 for downregulated).
Gene ontology and pathway analysis. Different Public reference databases including National Centre for Biotechnology Information (NCBI) non-redundant (nr), Swiss-Prot and UniProt Reference Clusters (UNIREF) were used to fix unigene identity through Basic Local Alignment Search Tool (BLAST) that has a sequence similarity index of up to an E value < 1.0E−5. Annocript program was employed for Unigene classification and assignment of GO terms to assembled transcripts. Functional categories were based on cellular components, biological processes, and molecular functions. Web Gene Ontology Annotation Plot (WEGO) tool was used for constructing GO plots representing sorted transcripts. Enrichment analysis of the top 20 DEGs was conducted. The Cluster of Orthologous Groups (COG) database was employed to achieve operational cataloging of unigenes and their role in different metabolic pathways was deduced through the KEGG database. DEGs of different comparison assemblies were represented employing numerous plots like HeatMap, Circos, and Map-Man (v3.51R2).
Quantitative real-time PCR (qRT-PCR) analysis. Total RNA was isolated from frozen leaf and root samples using Spectrum™ Plant Total RNA Kit™ (Sigma) according to the manufacturer's protocol and stored at − 80 °C. The quality and concentration of the isolated RNA were assessed by a NanoDrop spectrophotometer (Nano 200). The cDNA was synthesized from total RNA using (Takara Bio) as per manufacturer protocol and stored at − 20 °C. The coding sequence of selected DEGs was obtained from NCBI and gene-specific qRT-PCR primers were designed using IDT Primer designer and cross-checked by NCBI Primer-BLAST. The list of primers is presented in Supplementary Table S11. The qRT-PCR was performed in triplicates using the real-time PCR (Agilent Technologies, USA) detection system as described elsewhere 70 . The qRT-PCR of selected genes was carried out in those comparisons only (out of six total comparisons) in which the particular gene has significant expression in RNA-Seq analysis. The PCR program was set for 40 cycles. The maize Actin gene was used as an internal control to normalize gene expressions. Melting curves were analyzed and the relative fold change (log2 scale) in gene expression was calculated using the 2 −ΔΔCt method 71 .

Plant ethics statements.
All the methods used in the current study were carried out following relevant guidelines and regulations. All experimental protocols were approved by a named institutional and/or licensing committee.